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ABSTRACT 

The radiative cooling of optically thin gaseous regions and the formation of a two-phase 
medium and of cold gas clouds with a clumpy substructure is investigated. We demonstrate how 
dumpiness can emerge as a result of thermal instability. In optically thin clouds, the growth rate 
of small density perturbations is independent of their length scale as long as the perturbations 
can adjust to an isobaric state. However, the growth of a perturbation is limited by its transition 
from isobaric to isochoric cooling when the cooling time scale is reduced below the sound crossing 
time scale across its length scale. The temperature at which this transition occurs decreases with 
the length scale of the perturbation. Consequently small scale perturbations have the potential 
to reach higher amplitudes than large scale perturbations. When the amplitude becomes 
nonlinear, advection overtakes the pressure gradient in promoting the compression resulting in 
an accelerated growth of the disturbance. The critical temperature for transition depends on 
the initial amplitude. The fluctuations which can first reach nonlinearity before their isobaric 
to isochoric transition will determine the characteristic size and mass of the cold dense clumps 
which would emerge from the cooling of an initially nearly homogeneous region of gas. Thermal 
conduction is in general very efficient in erasing isobaric, small-scale fluctuations, suppressing a 
cooling instability. A weak, tangled magnetic field can however reduce the conductive heat flux 
enough for low-amplitude fluctuations to grow isobarically and become non-linear if their length 
scales are of order 10 -2 pc. If the amplitude of the initial perturbations is a decreasing function 
of the wavelength, the size of the emerging clumps will decrease with increasing magnetic field 
strength. Finally, we demonstrate how a 2-phase medium, with cold clumps being pressure 
confined in a diffuse hot residual background component, would be sustained if there is adequate 
heating to compensate the energy loss. 

Subject headings: globular clusters: general — instabilities — ISM: clouds, formation 
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1. Introduction 

The interstellar medium and cold gas clouds are characterized by a clumpy substructure and a turbulent 
velocity field (Larson 1981, Blitz 1993). As molecular clouds are the sites of star formation, their formation, 
internal structure and dynamics determines the rate of star formation and the properties of young stars, 
such as their mass function or binarity. The understanding of the origin of cold clouds and their internal 
substructure is therefore of fundamental importance for a consistent theory of star formation and galactic 
evolution. 

In the nearby clouds, the dispersion velocity inferred from molecular line width is often larger than 
the gas sound speed inferred from the line transition temperatures (Solomon et al 1987). MHD turbulence 
may be responsible for the stirring of these clouds (Arons & Max 1975). This conjecture is supported by 
the polarization maps and direct measurements of field strength in some star forming regions (Myers & 
Goodman 1988, Crutcher et al. 1993). Recent simulations of MHD turbulence, however, suggest that it 
dissipates rapidly (Gammie & Ostriker 1996, MacLow et al 1998, MacLow 1999, Ostrikcr ct al. 1999). One 
possible source of energy supply is winds and outflows from young stellar objects (Franco & Cox 1983, 
McKee 1989). But in regions where star formation is inactive, clumpy structure with velocity dispersion is 
also observed. Thus, the origin and energy supply of clumpy cloud structure remains an outstanding issue. 

On small scales, magnetic field pressure is important in regulating infall and collapse of protostellar 
clouds and the formation of low-mass stars (Mouschovias & Spitzer 1976, Nakano 1979, Shu 1993). For 
clouds with sub-critical masses, gravitational contraction is proceeded by ambipolar diffusion which for 
typical cloud densities operates on a timescale tb ~ 10 7 ~ 8 yr (Lizano & Shu 1989, Mouschovias 1991). 
In regions with intense star formation activities such as the central region of Orion, tb for individual 
dense clumps is comparable to the typical age of the young stellar objects. But, the spread in stellar ages 
(At* <~ 10 6 yrs) appears to be considerably shorter than tb (Carpenter et al. 1997, Hillenbrand, 1997). 
This coeval star formation history requires either a coordinated trigger mechanism for star formation within 
initially magnetically supported clumps or subcritical collapse, fragmentation and star formation of a larger 
molecular cloud region in which the magnetic field plays a weak role. 

A rapid and coordinated episode of star formation can also be inferred in globular clusters (Brown et 
al. 1991, 1995, Murray & Lin 1992, Lin & Murray 1992). In some metal deficient clusters such as M92, 
the total amount of heavy elements corresponds to the yield of a few supernovae. If star formation has 
proceeded over a duration At, comparable to the expected life span (<~ a few 10 6 yrs) of massive stars, 
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a significant metallicity spread would be expected, in contrast to the observations (e.g. Kraft 1979). At 
least in these systems, At* < tb and star formation may have proceeded through supercritical collapse. 
The dynamical timescale of most clusters at their half mass radius is Td ~ 10 6 yr. Any energy dissipation 
associated with the episode of star formation would imply an even longer dynamical timescale in the proto 
cluster cloud prior to that event. We infer that At* was comparable to or shorter than the dynamical 
timescale of the proto cluster clouds indicating a rapid fragmentation and star formation episode. 

In this paper, we focus on the rapid emergence of clumpy structure during the formation and collapse 
of a thermally unstable supercritical cloud. This process is relevant to the formation of stellar clusters as 
well as galaxies. We assume that the clouds condense out of a diffuse hot medium as a result of thermal 
instability. Large condensations with cooling timescales t c > Td are thermally stable because they can 
adjust through contraction such that their radiative losses may be compensated by the release of their 
gravitational energy. Runaway cooling of the gas through thermal instability however occurs in clouds with 
t c < Td- In order to form clumps within an initially almost homogeneous cloud, internal density fluctuations 
must grow rapidly on a timescale short compared to the mean dynamical timescale of the entire cloud. 
Small scale density fluctuations would begin to dominate if either the growth timescale or the limiting 
amplitude is a decreasing function of the perturbations' length scale. One possible fragmentation mechanism 
is gravitational collapse. The reduction in the cloud's temperature reduces its Jeans' mass, leading to the 
onset of gravitational instability and collapse. However, for a non rotating, cold, homogeneous gaseous 
region, gravitational instability alone cannot induce fragmentation because the growth rate is essentially 
independent of length scale such that the growth timescale for the density contrast is comparable to the 
dynamical timescale of the whole cloud (Hunter 1962). This has also been shown by numerical collapse 
simulations of initially gravitationally unstable perturbed gas clouds (e.g. Burkert & Bodenheimer 1993, 
1996, Burkert, Bate & Bodenheimer 1997). If the initial density perturbations So are linear (5 < 1), 
fragmentation is suppressed until the gas cloud has collapsed into either a disk or a dense filamentary 
substructure. 

We propose that clumpyness in clouds arises naturally from their formation through a cooling 
instability which acts on timescales that can be much shorter than the dynamical timescale of the cloud. In 
a pioneering paper, Field (1965) derived a criterion for a cooling gas to be unstable to the growth of thermal 
condensations. He showed that thermal instability can lead to the rapid growth of density perturbations 
from infinitesimal Sq to nonlinear amplitudes on a cooling timescale t c which for typical conditions in the 
interstellar medium is short compared to the dynamical timescale. If t c increases with decreasing density 
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any small density difference would induce a temperature difference between the cooler perturbed region and 
the warmer background. Across the interface between the two-phase medium, differential cooling leads to a 
pressure gradient which induces a gas flow from the lower-density background towards the higher-density 
perturbed region. The density enhancement in the cooler region further reduces its cooling timescale 
compared to that of the background where r c increases. A more detailed investigation of the growth of 
condensations in cooling regions has been presented by Schwarz et al. (1972) who included also the effects 
of ionization and recombination and by Balbus (1986) who examined the effect of magnetic fields. The 
classical model of the interstellar medium where heating balances cooling was presented by Field et al. 
(1969). A recent progress report on the theory of thermal instability is given by Balbus (1995). 

Although thermal instability proceeds faster than the collapse of the cloud, its growth rate is determined 
by the local cooling rate. During the initial linear evolution, variations in the initial over density (or under 
temperature) might lead only to a weak dependence of the growth timescale on the wavelength. In this 
paper we show however that there exist two important transitions which are very sensitively determined by 
the wavelengths of perturbations. 1) The growth of a perturbation is limited by its transition from isobaric 
to isochoric cooling, when the cooling time scale is reduced below the sound crossing time scale across the 
wavelength of the perturbation. This transition occurs at a lower temperature, with correspondingly larger 
over density, for perturbations with smaller wavelengths. 2) For those perturbation which can become 
nonlinear before the isobaric to isochoric transition, advection overtakes the pressure gradient in promoting 
the compression and growth of the perturbed region at an accelerated rate. The fluctuations which can first 
reach nonlinearity would dominate the growth of all perturbations with longer wavelengths and homogenize 
disturbances with smaller wavelengths. Thus, they determine the characteristic size and mass of the cold 
dense clumps which would emerge from the cooling of an initially nearly homogeneous cloud. Thermal 
conduction could in general erase these fluctuations, suppressing the instability. Weak, tangled magnetic 
fields would however be efficient enough in reducing the conductive flux, allowing the medium to break up 
into cold clumps on the characteristic length scale. 

We study the cooling and fragmentation of gas using simplified power-law cooling functions. Since we 
are primarily interested in supercritical clouds, we neglect the effect of magnetic fields. Note that even a 
weak magnetic field could have an important destabilizing influence in thermal instability (Loewenstein 
1990, Balbus 1995). In §2, we obtain approximate analytic solutions which describe the evolution of a linear 
density perturbation in the isobaric and nearly isochoric regime. We show that the growth of over density 
in a thermally unstable fluctuation is limited by a transition from isobaric to isochoric evolution and that 
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the limiting amplitude is a decreasing function of the length scale. We verify our analytic approximations 
with numerical, hydrodynamical calculations which are also used in §3 to study the transition into the 
non-linear regime. In §4 we investigate the cooling of interacting perturbations and determine the critical 
length scale of clumps that emerge through thermal instability. The importance of thermal conduction is 
investigated in §5. In §6 we discuss the affect of heating processes and the formation of a stable 2-phase 
medium. Finally, we summarize our results and discuss their implications in §7. 



2. The Initial Evolution of Thermal Instability 

The dynamical evolution of the gas is described by the hydrodynamical equations 

k=l 

k=l K fc=l K v 

where j=l,2,3 is the coordinate index, C v = R g /fi{T — 1) is the heat capacity, R g ,[i and T are the gas 
constant, mean molecular weight and adiabatic index, respectively. 

In the unperturbed state, the gas remains at rest (Uj = 0) and its density attains a constant value, p . 
The time dependent energy equation (3) gives 

F)T 

Cv ~dT = ~ PoA (4) 

where the cooling rate A = AoXjf . The power index is determined by the detailed atomic processes. Since 
we are primarily interested in the physical evolution of thermal instability, we adopt a simple constant [3 
prescription. The cooling would be thermally unstable (with r c = To/dTo/dt as an increasing function of 
To) in the isochoric region if (3 < 1 and in the isobaric region if /3 < 2. In the absence of external heating, 
the unperturbed gas temperature T can be expressed as a function of the dimensionless time variable 
t = t/r c (0) such that 

T (i)=To(0)(l-(l-/3)r)^ (5) 

where T o (0) and r c (0) = C v / / ooAoT (0)' 3 ~ 1 are the initial (at t = 0) temperature and cooling timescale, 
respectively. 
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2.1. The Perturbed Quantities 

The evolution of the perturbed density (pi = p — po), temperature (Ti = T — To), and velocities (Uj) 
are derived from the linearization of the equations (1) to (3): 

lei = _Y?Ei (6) 

dt po f-f dx, ' v ; 

dUj _ R g Tp d (T x | gi 
at p po 

and 



(7) 



otT f-f ctej t c \p T J 



where 



is the characteristic cooling timescale at the instant of time t. 

Since the perturbation equations are linear in Xj, we adopt a local approximation in which the 
positional dependence of all the perturbed quantities is proportional to exp(ikjXj) where kj is the wave 
number in the jth direction. Substituting a dimensionless velocity variable Vj = ikjT c (0)Uj, the perturbed 
equations reduce to 

dPl =-±V, (9) 



drpo J=1 



dVj 



(10) 



8r 

where ^ = |f + ^ is the perturbed pressure, Kj = T c (0)kjy/ R g To (0) / p is the ratio of the initial cooling 
to sound crossing timescale over a characteristic wavelength 2ir/kj, and 

For a perfect gas, the unperturbed pressure Po = RgPoTo/p decreases at the same rate everywhere. We find 
from Eqs ([)]) and (|ll|) that the amplitude of the perturbed pressure is 

3 



±*L = ±( * + -rTv, L 

drP Or \p T ) ~ ^ 3 1 - (1 - 



(2-/3)^-(l-/3)^]. -12) 
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2.2. The initially isochoric regime with K < 1 

For computational simplicity, we now consider a 1-D limit treatment in which the initial (at r = 0) 
amplitude of pi equals to a finite value p a with that of V± and Pi equal to zero. These conditions correspond 
to an initially almost homogeneous, hot region of gas in pressure equilibrium. To third order in r the eqs 
(0), (0), and ( |l2|) give the following approximate solution 

^^ P -(l + ^(2-P)rA (13) 
Pa Pa V 6 / 

V „ (P- 2)if 2 Pa , 3t 2 + 2/3t3 ^ (u) 
6 po 

p Ktf _^ + (1 _ ffl ^ + (|_^_™!),j) (15) 

Figure 1 compares this solution with a numerical integration of the complete non-linear hydrodynamical 
equations (1) to (3) for K=l and K=0.5. We use a 1-dimensional version of the second-order Eulerian 
hydro code which is described in Burkert and Bodcnheimer (1993). The agreement between the numerical 
results (solid lines) and analytical solution (dots) is excellent, even for large values of r ~ 1 where the basis 
of the analytic approximation is no longer valid. 

Due to slightly more efficient cooling within the density perturbation a small pressure gradient builds 
up. In an attempt to maintain pressure balance, the slightly warmer gas in the low-density regions 
continually compresses the more dense and cooler parts. Consequently the over-density in the perturbed 
region increases as the gas cools. Figure 1 and the equations (13) to (15) show that for small values of 
K < 1 the growth rate of the density enhancement depends on the size of the perturbation and increases 
with increasing values of K or decreasing wavelength. However, due to its long sound crossing timescale, 
the perturbation cannot be compressed significantly while cooling; it cools almost isochorically (Parker 
1953). After one cooling timescale the gas temperature has reached its minimum value with the density 
enhancement still in the linear regime. Now, the pressure gradient reverses, erasing the fluctuation. 

2.3. The initially isobaric regime: K > 1 

The solid lines in figure 2 show a numerical calculation of the evolution of a density perturbation with 
K = 200, (3 = and T = 5/3. For K >> 1, perturbations can react quickly on any pressure gradients due 
to the short sound crossing timescale, relative to the cooling timescale. The simulations indicate a solution 
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which consists of a fast oscillatory part and a slowly growing part. Linearizing the slowly growing part, we 
find from the equations (9) to (12) the following approximate solution: 



" r ^^,,1 (16) 



""^fl^^TTi + '"')• (17) 



po po V — 2 + /3 ) \ r 
-2 Pa ( 

~po~ ^^r-2 + /3 

and 

ft . t^J™ * ) (1 + (1 " & ~ ^ 

The characteristic frequency is determined by a cubic dispersion relation 

lu 3 + i(l - I3)lj 2 - K 2 Tlu - i(2 - P)K 2 = 0. (19) 

A similar relation was discussed by Balbus (1995). Note that the linearized solution is also valid for r > l/u> 
as long as t << 1. As K >> 1 the two dominant real roots {lo » ±y/TK) of the dispersion eq([l9|) yield 
oscillatory parts in px/po, V, and P\/Pq. The upper panels of figure 2 show that for r < 0.1 the analytical 
approximation is in good agreement with the numerical integration of the nonlinear hydrodynamical 
equations. 

For all length scales, the ratio of sound propagation to cooling timescale 



Q = T c (t)kyjR g T (t)/n = K(l - (1 - /3)t)3=$ (20) 

decreases during the subsequent evolution. Provided Q » 1, cq(|Tg|) implies that the magnitude of Pi/Pq 
is much smaller than both V and pi/po- That is, the fluctuation reacts isobaric. Adopting P\/Pq ~ and 
neglecting the oscillatory term, the equations (9) and (12) can be combined to 



-V=±(H) = * (21) 

dr\ P J r(l-(l-/3)r)p 1 1 



with the solution 



and 



Pl = P±(l-(l-P) T )-TX=W (22) 
Po Po 



V = - ( 2—Ae±{i-{i-p) T )-T^mT-\ (23) 
r po 
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In contrast to fluctuations with K < 1 the evolution of isobaric fluctuations is independent of K. For 
(1 — >> t > solutions for pi/po in Eqs ( |l6| ) and ( p2f ) are in agreement to first order in r. The 
lower panels of figure 2 show that within a cooling time both pi/po and V are amplified to very large 
values. The agreement between the numerical calculation and the analytical prediction (equation 22 and 
23) is excellent. The opposite signs of pi/po and V confirm that mass is being pushed into the cool dense 
regions. For r « 1, the numerically derived density enhancement falls below the predicted values as the 
fluctuation becomes isochoric and contributions from the perturbed pressure cannot be neglected anymore. 



2.4. Transition to Isochoric Evolution and the Emergence of Small Scale Perturbations 

Figure 3 shows the density evolution of initially isobaric fluctuations with different ratios of cooling 
to sound crossing times K as determined from the numerical calculations. The initially isobaric growth of 
the density fluctuations is independent of wavelength and K and in excellent agreement with equation (22) 
(dashed curve). The perturbations transform however to the isochoric solution for the epoch after Q has 
declined below unity. Thereafter, gas in the perturbed region cools off faster than it can adjust to a pressure 
equilibrium with the surrounding region. Subsequently, the over density of the perturbed region is slowly 
modified by the inertial motion Vtrans of the gas at the time of the transition and its growth stalls. In 
Figure 3 the transition into the isochoric regime is indicated by the overdensity falling below the expected 
value shown by the dashed thick line. 

Although the growth of the perturbed quantities does not explicitly depend on the wavelength fc, the 
critical transition time when Q ~ I 

2/3-2 

1 - 
i - p 

is a function of K (and k). At this transition point, the over density in the perturbed region is 

^i = «. (25) 
Po Po 

and the velocity is 



v trans =- 2 -^A p -i K &^^i (26) 

r po 

In the isochoric regime the amplitude of the perturbed density increases as 
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Pi 



Ptrans 




(27) 



which increases much less steeply than the isobaric fluctuations (equation 22). 

For thermally unstable clouds, (3 < 1 such that ptrans/po is an increasing function of K or a decreasing 
function of the wavelength (A) of the perturbations. Despite the independence of the rate of change of pi/po 
on A, equation 24 shows that for (3 < 1 the short length scale disturbances undergo isobaric to isochoric 
transition at a later time and therefore acquire a greater limiting amplitude than the long length scale 
disturbances. Thus, the short length scale disturbances would emerge to dominate the structure of the cloud 
unless the initial perturbation amplitude p a /po increases with K more rapidly than lf( 2 ^~ 4 )/( 3 ~ 2 ' 3 ) r . This 
evolution is physically equivalent to the fragmentation process in which the contrast between the enhanced 
density in a disturbance and the average cloud density becomes most pronounced on the smallest scales. 



In Fig. 3 fluctuations with very large values of K show yet another evolution: for later times the 
overdensity rises faster than predicted by equation (22). These fluctuations become nonlinear with 
Pi / Po > 1 before the transition into the isochoric regime. The critical value of K for this evolution can be 
estimated from equation (25) assuming ptrans/po = 1 ; 



For K > K cr it the analytical approximations discussed previously are not valid anymore and we have 
to investigate the evolution numerically, solving the complete non-linear hydrodynamical equations. The 
simulations shown in figure 3 assumed /3 = 0, T = 5/3 and p a /po = 10 -3 . For these values the simple 
approximation (28) predicts K cr n = 5600 which is roughly in agreement with the numerical results where 
the transition into the nonlinear regime occurs more smoothly between K=1000 and K=5000. Equation 
(28) somewhat overestimates K crit because nonlinear effects actually become important earlier, when the 
overdensity is in the range 0.1 < pi/po < 1- 

Figure 4 shows the structure and evolution of a non-linear fluctuation. During the early isobaric 
evolution the pressure gradient (lower right panel) is negligible. A small pressure gradient builds up in 



3. The Transition into the Nonlinear Regime 




) 



(2f3-3)r 
4.-2/3 



(28) 
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the nonlinear regime, where the profiles cannot be approximated anymore by sinusoidal functions but 
instead become strongly peaked towards the center. Non-linear fluctuations grow fast with the density and 
temperature reaching their maximum and minimum values, respectively, at a time T cr n < 1 which is shorter 
than a cooling time. Due to the fast growth in the nonlinear regime, t ct h is roughly given by the time when 
Pi/ Pa = 1- From Eq (^), we find 

1-/3 r Vpo, 



T cH t = ^— a 1-f-l I- 



At t = T cr it, the dimensionless velocity (see Eq. |23|) 



V crit = V(r crit ) = ^— ? f ^ ) " (30) 



r VPo, 

is much larger than unity for perturbations with small initial amplitudes such that contributions due to 
nonlinear advection (such as Ujdp/dxj, UjdUj/dxj, and UjdT/dxj) would exceed the linear contributions 
contained in the perturbed equations ([||8|) before the over density p\ has become comparable to po (see 
above). Advection generally enhances the effect of compression and promotes the growth of density contrast 
at an accelerated rate. 

Although the time of maximum compression for a fluctuation with K > K cr n does not depend 
explicitly on the length scale, it is determined by the initial amplitude p a / Po of the perturbation which may 
be a function of the wavelength. For an initial power-law perturbation in which p a /po — Ao(k/ko) ri , 

1 / / 1 fc \ 2-/J 

'AS-) . (31) 



If the amplitude of the initial perturbation is an increasing function of the wavelength (which 
corresponds to a negative 77) , t ct h would be an increasing function of k in the thermally unstable region 
with p < 1. In this case, nonlinearity would be first reached on the largest length scale with K > K cr n. If, 
however, the amplitude of the initial perturbation is a decreasing function of the wavelength (i.e. rj > 0), 
Tcrit would be a decreasing function of k and nonlinearity would be reached on the smallest scale first. Note 
that for 1 < j3 < 2, the dependence of r cr i t on 77 and k is reversed. 

In Figure 3, a temperature independent cooling function has been used. In order to determine the 
dependence on the specific form of the cooling function, additional simulations have been performed, 
adopting a more realistic cooling function (Dalgarno & McCray 1972) which assumes solar element 
abundance and collisional equilibrium ionization. Note that for temperatures T > 10 4 K the cooling rate 
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is several orders of magnitudes larger than for T < 10 4 K, defining two different temperature regimes with 
very different cooling timcscales. The simulations show that the previous results remain valid for each of 
these temperature regimes. Starting in the low-temperature regime, a fluctuation will become non-linear 
for K > K cr i4 and cool down to the lowest allowed temperature. The same is true for fluctuations that start 
in the high-temperature regime. Non-linear fluctuations in this regime do however stop cooling efficiently 
at T w 10 K, leading to high-density clumps with such a temperature. 

4. Interacting Fluctuations and the Emergence of Substructure with a Critical Wave Length 

Up to now we have investigated the evolution of isolated fluctuations. In reality however a cooling 
gaseous region consists of a superposition of fluctuations with different wavelengths and amplitudes. As 
we indicated in the previous section, the outcome of the thermal instability may be determined by the 
wavelength dependence in the initial amplitude of the perturbations. 

In order to illustrate various competing effects such as isobaric to isochoric transition and the onset 
of nonlinear growth, we present in Figure 5 a series of models with (3 = and T = 5/3, where the initial 
density distribution consists of the superposition of two fluctuations with ratios of wavelengths A1/A2 = 20 
and amplitude ratios p a ,\l Pa.i = 2 which corresponds to rj = —0.23. Four values of Ai were chosen and 
they correspond to Xi = l,10,100 and 1000, respectively. The K2 values for the smaller perturbation are 
always a factor 20 larger. Since the initial overdensity p a ,i/po — 0.01, the critical value of K for nonlinear 
evolution is K cr n = 316, according to equation (p8|). We show the density distribution after 1 r c (0). In the 
upper left panels of Fig. 5, the fluctuations have values of K\ = 1 and K2 = 20 which are small compared 
to K cr i t . Their growth therefore stalls due to transition into the isochoric regime and the overdensity after 
a cooling time is still linear. The smaller fluctuation dominates at the end because its isochoric transition 
occurs later and at a higher overdensity than for the larger perturbation. In the upper right panels with 
K\ = 10 and K2 = 200 the smaller perturbation is again dominating after r = 1 although, now, the density 
distribution is also affected by the underlying larger perturbation. In both cases the density within the 
density peaks does not decrease much with respect to its initial value. The situation is different in the lower 
left panel with K\ = 100 and K% = 2000. Here the smaller perturbation has become nonlinear, generating 
small dense clumps which stand out against the larger perturbation. Up to now, the smaller perturbation 
was always dominating the density distribution after a cooling time. The situation is however different in 
the lower right panel where K\ > K cr n. Now, the larger perturbation becomes nonlinear and advection 
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drives all the gas and its small scale fluctuations into one very dense, cold clump that is embedded in a hot 
diffuse environment, erasing smaller scale fluctuations. 

The dependence of structure formation on the initial power-law perturbation index r\ is illustrated 
in figure 6 which shows the initial and final density distribution of two interacting perturbations with 
A1/A2 = 10, K\ = 2 x 10 4 , p a ,i/Pa,2 = 10 (77 = —1) in the left panels and p a ,i/p a ,2 = 0.1 (r/ = 1) in the right 
panels. As expected, in the case of 77 = —1 the larger perturbation becomes nonlinear first, leading to one 
massive density peak after a cooling time. For 77 = 1, the small length scale perturbations begin to dominate 
after a cooling time, breaking the region up into dense clumps on the smallest scale. More specifically, if 
the amplitude of the initial perturbation increases with increasing wavelength (77 < 0), clumps will form 
with length scales A w \ cr it- Otherwise, the sizes of the fastest growing perturbations will be A w X K . In 
this case, the clump sizes should decrease with increasing magnetic field strength. 



During the growth of linear density perturbations in the isobaric regime the resulting temperature 
gradient will induce conductive heating of the fluctuations. 



Several studies (e.g. McKee & Begelman 1990, Ferrara & Shchckinov 1993) have demonstrated that 
thermal conduction could stabilize and even erase a density perturbation if its scale is smaller than the 
Field length (Field 1965) 



where k is the thermal conduction coefficient. 

In order to include the effect of thermal conduction, the term V (kVT) / pC v has to be added to the 
right-hand side of equation (3). The linearized pressure equation (12) is then 



5. The importance of thermal conduction 



5.1. Thermal conduction in the absence of magnetic fields 




(32) 



d_F\ 
drP 



3 



1 



-(!-/*)£) 
V Po PqJ 




l-(l-/3)r 



fc 2 T c (0). 



(33) 
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In the isobaric regime with P\/Pq ~ and Ti/T w —pi/po thermal conduction will become important 

if 

In the early stages of cooling (r << 1) fluctuations will therefore be erased by thermal conduction if their 
wavelenghts are 

2ir 

Figure 7 shows the evolution of an initially isobaric, 1-dimensional fluctuation with a ratio of cooling- 
to sound crossing time K=200 and (3 = 0. The 1-dimensional, non-linear hydrodynamical equations (1) 
- (3) are solved numerically, including thermal conduction. The solid line shows the evolution of the 
density contrast pi/po as predicted by the analytical model (equation 22) which is in excellent agreement 
with the numerical result (filled points) for A > A K . For A = A K (upper dashed line) conductive heating 
is non-negligible anymore and the fluctuation grows less fast. For A < A K the growth of fluctuations is 
suppressed by thermal conduction. 

In summary, thermal conduction can play a significant role in regulating the break-up of a radiatively 
cooling gaseous medium. Small scale substructure can only emerge in a limited wavelength regime which is 
given by 



A K < A < X cr it (36) 

The growth of perturbations is completely suppressed if \ cr it < A K for all perturbations. 

The dotted lines in figure 8a show \ crit for fluctuations with initial overdensities 
\og(p a /po) = —3,-2,-1, the solid line shows A K . A cooling function assuming collisional equilibrium 
ionization (Spitzer 1978, Dalgarno & McCray 1972) and a realistic conduction coefficient (Ferrara & 
Shchckinov 1993) has been adopted. The dashed line shows the mean free path (Cowie & McKee 1977) 



(37) 
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for electron energy exchange. Note that for a given temperature the ratios X K /X cr i t / X e are independent of 
pressure. 

The classical thermal conductivity is based on the assumption that A e is short compared to the 
temperature scale height hx ~ \/{T\/Tq) > X cr it- Otherwise, the heat flux q is saturated (Cowie & McKee 
1977) and no longer equal to q = — kVT. Indeed, figure 8a shows that linear fluctuations in astrophysical 
plasmas will in general lie in the non-saturated regime (X cr it > X e ) for initial amplitudes p a /po > 0.001. 

However, we also find that in general X K > X cr i t for perturbations with amplitudes log(p a /p ) < —2. 
This implies that a cooling instability, resulting from linear density perturbations will in general be 
suppressed by thermal conduction. 



5.2. Thermal conduction, including magnetic fields 

The interstellar medium is in general penetrated by magnetic fields. In most situations the electron 
mean free path X e is large compared to the length scale at which the resistive destruction of the magnetic 
field is significant. A tangled magnetic field can then develop, concentrated on scales Is which are smaller 
than A e (Chandran & Cowley 1998). When the gyroradius 

of thermal electrons with typical velocities vt = (kT/rrie) 1 / 2 is much smaller than Ib or A e the magnetic 
field controls the motion of individual electrons. This condition is satisfied in many astrophysical plasmas 
even if the magnetic field is too weak to be hydrodynamically important. 

If a is small compared to the length scale A of a fluctuation, heat is conducted according to the classcal 
thermal conduction equation (Spitzer 1962), however with a thermal conductivity Kb which is reduced from 
the classical Spitzer value k as a result of the tangled magnetic field by (Chandran & Cowley 1998) 

ln(< B /a) 

If we normalize B to its value for magnetic-to-thermal energy equipartition B T = {2AirpR g T) 1 / 2 we 

find 
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For weak magnetic fields (B w 0.01-Bt) and length scales of order the electron mean free path 
equation 40 leads to ln(Zs/a) > 10, by this reducing Kg by two orders of magnitudes and the length scale 
of thermal conduction by one order of magnitude. As an example, the shaded area in figure 8b shows the 
wavelength regime A K < A < \ cr it where linear fluctuations with p a /po = 10~ 2 could grow as a result of 
cooling, assuming a pressure of P/ks = 10 3 K cm~ 3 and Ib — 0.1A e . The presence of a weak magnetic field 
can suppress thermal conduction efficiently, allowing small scale structure with wavelengths A w 10~ 2 pc to 
emerge as a result of cooling. 



6. The importance of heating and the emergence of a stable two-phase medium. 

The calculations in sections 3 and 4 showed that cooling gas clouds with small initial perturbations 
break up on a critical wave length X cr it below which over density first becomes nonlinear. If the initial 
amplitude is a decreasing function of A, the clouds would break up on the smallest length where the local 
radiative cooling law remains valid and fluctuations are not destroyed by conduction. But, if the initial 
amplitude is an increasing function of the wavelength, 



2ttt c (0) R g T (0) 
Acrit = \ • (41) 

small perturbations on scales A < \ crit are erased when the gas accumulates in the center of fluctuations 
with A = Xcrit- Perturbations with A > \ cr it do not become nonlinear but they break up into substructures 

with A = Xcrit- 

After a cooling time the dense, cold, non-linear perturbations are embedded in a warmer, diffuse 
environment (see Fig. 9). However, the example in figure 9 shows that the cooling timescale of the 
inter-clump gas remains short compared to the initial cooling timescale. This gas therefore cools to a ground 
state temperature T min shortly after r = 1. Subsequently the reversed pressure gradient would remove the 
fluctuations unless they are gravitationally bound. 

In order to maintain a stable two-phase medium a heating term must be included (Field et al. 1969). 
Here we assume a power law dependence of the heating rate 
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T h = Top 1 



(42) 



where T and 7 are constants. In general, the size of the whole cooling region is large compared to X crit 
such that it cannot establish pressure equilibrium with the surrounding confining medium during a cooling 
timescale. In this case, the region cools isochorically and breaks up into substructures on scales of A crit 
before establishing pressure equilibrium with the environment. If the average gas density po is smaller than 
a critical value 



heating would dominate everywhere and the region would adjust to a thermal equilibrium state where 
heating is balanced against cooling. If po > pr cooling dominates and the density fluctuations would 
grow and become non-linear as discussed in the previous sections. Eventually after a cooling time, 
the region would break up into cold high-density condensations which are separated by warm gas with 
densities p m in « Po- If Pmin > Pr this interclump medium would cool as shown in figure 9 and the 
density fluctuations would be erased. Figure 10 shows a situation with p m i n < pr- Heating dominates 
in the interclump region where the gas temperature and gas pressure rise, until pressure equilibrium is 
established. A stable 2-phase medium has formed with cold clouds of minimum temperature embedded in a 
hot interclump medium with a temperature that is determined by the balance of cooling and heating. 



The discussions in this paper focussed on the emergence of small scale perturbations. We have assumed 
the pre-existence of small initial perturbations which is a reasonable assumption for dynamically evolving 
systems like the interstellar medium in galaxies or in galactic clusters. We have limited our analysis to the 
optically thin regime such that radiation transfer is solely due to optically thin local radiative processes. 
This approximation is appropriate for the collapse of supercritical clouds where the effect of a magnetic 
field is dominated by thermal processes. Such a situation may be particularly relevant for the formation of 
stellar clusters and first generation stars in galaxies. Provided that the density of the progenitor clouds is 
relatively small, the local cooling approximation is adequate. We also neglected the interaction and merging 
of clumps. These processes become important for the subsequent evolution and they will be considered in 



1 




(43) 



7. 



Discussions 
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subsequent papers. 

In the context of our approximations, we have shown that thermal instability can lead to the breakup 
of large clouds into cold, dense clumps with a characteristic length scale which is given by X cr it in eq. (^TJ) 
or by the smallest unstable wavelength that is not erased by thermal conduction, depending on whether 
the amplitude of the initial perturbation is an increasing or decreasing function of wavelength. For linear 
perturbations with overdensities p a /po ~ 0.01 the critical wavelength lies in the regime of 10 -3 pc to 10 _1 
pc, depending on the initial temperature. The emergence of small scale dense subcondensations is equivalent 
to fragmentation. As in a thermally unstable region the cooling timescale is shorter than the dynamical 
timescale, gravity has no time to play an important role during this fragmentation process. \ cr it may be 
either smaller or bigger than the Jeans' length. In the latter case gravity becomes important eventually. In 
general however, thermally induced fragmentation of clouds with small initial density fluctuations proceeds 
the onset of gravitational instability of their individual clumps. 

In our analyses, we adopted an idealized power-law cooling function. In reality, the cooling efficiency 
would terminate when the main cooling agents reach their ground state or establish an equilibrium with 
some external heating source. The latter is necessary for the clouds to attain a two-phase medium. 
Interaction between these two phases may determine the pressure, density and infall rate of the cloud 
complex as well as the dynamical evolution and size distribution of cloudlets and sub condensations. The 
analysis of this interaction will be presented elsewhere. 
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Fig. 1. — The solid lines show the numerical results of the evolution of the overdensity (left panel) and 
dimcnsionless velocity (right panel) for almost isochoric perturbations with K=l (upper curves) and K=0.5 
(lower curves). A cooling law with (3 — is adopted. The dots represent the analytical solution (equations 
13 and 14). 
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Fig. 2. — The growth of an initially isobaric density perturbation with K=200 is shown, adopting (3 = and 
r = 5/3. The numerical results (solid lines) are compared with the analytical solutions (dashed curves). The 
upper two panels show the early evolution which is compared with the predictions of the linearized equations 
(16) and (17) assuming uj = \/TK. The lower panels show the evolution till r = 1 which is compared with 
the results of the nonlinear equations (22) and (23). 
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Fig. 3. — The growth of initially isobaric fluctuations is shown. The solid curves, labeled by the adopted 
initial K-valucs of the fluctuations show the results of numerical calculations for (3 = 0. The thick dashed 
line shows the evolution of a linear isobaric perturbation as predicted by equation (22). Fluctuations with 
small values of K become isochoric and their density increase stalls. Above a critical value of K fluctuations 
become nonlinear during their isobaric growth. These fluctuations evolve faster than predicted by the linear 
approximation and achieve very large density contrasts on a timescale which is independent of K and shorter 
than a cooling timescale. 
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Fig. 4. — The evolution of overdensity pi/po, dimensionless velocity V, temperature Ti/T and pressure 
Pi/Po of a nonlinear fluctuation with K=1000, assuming (3 — 0. Cooling is assumed to stop after the 
temperature has decreased by 4 orders of magnitude. The evolutionary state of the system is shown at the 
time when the overdensity is 0.1, 0.5, 1 and 5 times the initial density. 
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Fig. 5. — The figures show the overdensity and the temperature distribution after a cooling time of regions 
with two interacting fluctuations with wavelength ratios Ai/Ao = 20 and amplitude ratios p a ,i/pa,2 = 0.5, 
adopting different values of K. The x-coordinate is normalized to the wavelength of a perturbation with 
K=l. 
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Fig. 6. — The left upper and lower panels show, respectively, the initial (r = 0) and final (r = 1) density 
distribution of two nonlinear interacting density perturbations with Ai/Ao = 10 and r\ = —1. The initial and 
final density distributions for interacting nonlinear fluctuations with r] = 1 are shown, respectively, in the 
upper and lower right panels. 
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Fig. 7. — The effect of thermal conduction on the growth of a perturbation with p a /po = 0.01 and K = 200 
is shown, adopting a cooling function with = 0. The solid thick line shows the theoretical prediction 
and the filled points show the numerical result for A K = 0. The lower dashed lines show the evolution of 
fluctuations with A = A K , A = 0.75A K , A = 0.5A K , and A = 0.25A K , respectively. 
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Fig. 8. — Fig. 8a shows several interesting lengthscales as function of temperature, adopting a gas 
pressure of P/ks = nT = 10 3 Kcm~ 3 . Solid line: A K . Dashed lines: X cr it for density fluctuations with 
log(p a //°o) = 3, respectively. Dashed line: A e . The shaded area in Fig. 8b shows the wavelength 

regime X K < A < \ crit where density perturbations with log(p a /po) = ^2 would be able to grow and become 
non-linear, adopting a gas pressure of P/ks = 10 3 Kcm~ 3 , a magnetic field of B = O.OIBt and a tangled 
magnetic field length scale of Ib = 0.1A e . 
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Fig. 9. — Evolution of a nonlinear perturbation with K=1000 and (3 = 0. The upper panels show the evolution 
of the maximum density and temperature (solid line) and of the minimum density and temperature (dashed 
lines) for 2 cooling timescales. Shortly after one cooling timescale the low-density regions cool down to the 
minimum temperature too and the reversed pressure gradient erases the density contrast. The lower panels 
show the evolution of the density and temperature during the short epoch atrwl when a 2-phase medium 
has formed. 
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Fig. 10. — The evolution of the maximum and minimum density (left panel) and maximum and minimum 
temperature (right panel) of a non-linear perturbation, including a cooling term with (3 = and a heating 
term which dominates for p < pr = O.lpo, where po is the initial average density. The density perturbation 
grows and becomes non-linear within a cooling timescale. As gas is pushed into the high-density regions, the 
density in the inter-clump region decreases below the ciritical value where heating begins to dominate and 
the inter-clump gas heats up. A stable two-phase medium forms. 



